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The thermal convection equations for a thin layer of fluid are solved numerically as an initial value 
problem. The calculations include only those nonlinear terms which have the form of an interaction of a 
fluctuation in the velocity and temperature with the mean temperature field. In the present calculations, 
the velocity and temperature fluctuations have one horizontal wave number, and satisfy free boundary 
conditions on two conducting horizontal surfaces. 

The computed steady state velocity and temperature amplitudes show many of the observed qualitative 
features. In particular, the experimentally observed boundary layering of the mean temperature field is 
correctly reproduced, and, at large Rayleigh number, the total heat transport is found to be proportional 
to the cube root of the Rayleigh number, provided the fluctuating temperature and velocity amplitudes have 
that horizontal wave number which maximizes the total heat transport. However, the heat transport found 
here for free boundaries is about three times the experimental value for rigid boundaries. The mean tem- 
perature gradient can become negative near the boundaries for large Rayleigh numbers and large horizontal 
scale motions. 

The linear stability of the system is also investigated, and it is concluded that the stable solutions for all 
Rayleigh numbers investigated (R<10®) have horizontal wave numbers which very nearly maximize the 
total heat transport. The stability study also indicates regions in which two or more horizontal wave numbers 
are required to support convection. ^ ^ ^ ^ 
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1. Introduction 

This paper describes the results of a numerical in- 
vestigation of the thermal convection equations for a 
thin layer of fluid confined between two plates on which 
free boundary conditions are employed. Our theoretical 
procedure is to include only those nonlinear terms which 
describe the interaction of the mean temperature with 
velocity and temperature fluctuations. That is, those 
terms responsible for eddy viscosity and eddy conduc- 
tivity effects on the turbulence itself are omitted. The 
above eddy terms (hereafter referred to as fluctuating 
self-interactions) are discarded in a physically consistent 
manner, so that no unrealistic behavior results from 

their omission. 

The motivation for this research is to discover quanti- 
tatively to what extent the turbulent convection prob- 
lem can be comprehended without the fluctuating self- 
interactions. The system of equations obtained by 
deleting these terms corresponds to closing the hierarchy 
of moment equations at the first nontrivial level by 
discarding third order cumulants. The resulting system 

of equations is complete and involves no empirical 
parameters. Moreover, the gross energetics of the flow 
are preserved. 

The method of numerical solution consists in inte- 
grating the Fourier amplitudes of the velocity and tem- 
perature fields forward in time until the steady state is 
achieved. This procedure has the advantage of assuring 


the stability of the final state provided a sufficient range 
of initial data is sampled. The present calculations, 
carried out on an IBM 7090 computer, contain one 
horizontal wave number and enough vertical wave 
numbers to ensure the elimination of truncation errors. 
In the numerical analysis, we have set the Prandtl 
number equal to unity. However, inspection of the 
equations which omit the fluctuating self -interactions 
shows that the heat transport is not a function of 
Prandtl number, if the system is steady. 1 

The calculated velocity and temperature fields show 
many of the qualitative features of the experimentally 
determined fields. In particular, at large Rayleigh num- 
ber, R, the total heat transport is found to be propor- 
tional to R*, provided the velocity and temperature 
fluctuation fields have that horizontal wave number 
which maximizes the heat transport. However, the heat 
transport found here for free-boundary conditions is 
three times the experimental value for rigid -boundary 
conditions. Preliminary numerical studies of the rigid- 
boundary problem indicate that for large Rayleigh 
numbers (R^ 10 6 ) the heat transport is about a factor 
of 2.3 smaller than that for free-boundaries and there- 
fore approximately 50 per cent higher than the experi- 
ment. Thus, it appears that the boundary conditions 
are quite important in producing the experimental heat 

1 For a discussion of the Prandtl number dependence of the heat 
transport for the complete system, see JCraichnan (1962a), 
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transport. The system has two additional failings: it 
turns out that the fluctuating amplitudes are steady in 
time and the horizontal plan form of the motions is 
indeterminate. 

The mean temperature gradient at low R closely 
resembles the experimental temperature profiles. At 
large Rayleigh numbers (R~10 6 ), the gross features of 
the temperature profiles are correctly predicted by our 
system. The computed mean temperature gradients are 
large in a thin layer adjacent to the boundary and are 
quite small in the body of the fluid. The gradients near 
the boundary can become negative for motions of large 
horizontal scale, but remain positive for motions of a 
sufficiently small horizontal scale. 

The stability of the steady state solutions against 
infinitesimal perturbation at other horizontal wave num- 
bers is also investigated and the regions of instability 
are delineated. These results closely parallel perturba- 
tion results at low Rayleigh number and support the 
idea that the most stable solution is near the one trans- 
porting the most heat flux (Malkus and Veronis, 1958). 
A finite amplitude stability study, and the associated 
development of a several-horizontal-wave-number sys- 
tem to steady state will be the topic of a future 
investigation. 

The idea that the important features of the convec- 
tion problem are contained in the system which omits 
the fluctuating self-interactions is implicit in the theory 
of convection advanced by Malkus (1954). In the origi- 
nal formulation of his theory, Malkus sought a maxi- 
mum for the heat transport subject to the constraint 
that the temperature gradient be positive, and that 
there be a smallest scale of motion participating in the 
advective heat transport. The smallest scale is supposed 
to be determined by the requirement that it be margin- 
ally stable in the presence of the mean temperature 
gradient occurring in the fluid. The smallest scale so 
determined furnished a cut-off in the cosine representa- 
tion of the mean temperature gradient. The assumption 
that the heat transport was maximum under the above 
constraints then led to an explicit form for the tempera- 
ture gradient. 

A more recent formulation of the Malkus theory by 
Spiegel (1962) replaces the cosine representation of the 
temperature gradient by an expansion in terms of the 
set of eigenfunctions, which are marginally stable on the 
mean temperature gradient. This version of the Malkus 
theory is exactly equivalent to the system considered 
here, provided the horizontal scale of the motions is such 
that the mean temperature gradient is everywhere 
positive. In this sense, our numerical results contain, 
as a special case, the exact solutions to the Malkus 
theory for one horizontal wave number and free 
boundaries. 

In this connection, a comparison of our computed 
temperature gradients with the predictions of the 


Malkus theory is relevant. In making this comparison, 
we must keep in mind that the system considered here is 
explicitly confined to only one horizontal wave number, 
whereas Malkus makes no explicit references to the 
nature of the horizontal-wave-number spectrum. We do 
not confirm the z~ 2 law for the gradient outside the 
boundary layer as predicted by Malkus, nor do we find 
a sharp cut off in the cosine spectrum of the temperature 
gradients. 

At low Rayleigh numbers (R<2000) our numerical 
results are in agreement with the calculations of Malkus 
and Veronis (1958) and Kuo (1961), who have obtained 
perturbation solutions to the convection equations. A 
procedure similar to ours has been used by Saltzman 
(1961) in studying the complete convection equations 
for R<6000. Our approach differs from his in that we 
are able to allow very many vertical modes to be excited, 
whereas his results are limited to one vertical mode. Our 
results indicate that it is essential to allow many more 
vertical modes than horizontal modes, if large trunca- 
tion errors are to be avoided. Thus, at R=4000, 10 
vertical modes must be included to obtain realistic 
temperature profiles. 

2. Theory 

a) The equation of motion and boundary conditions . We 
consider a thin layer of fluid confined between two 
infinitely conducting plates located at z = 0 and z=dr 
The lower plate is maintained at zero degrees, and the 
top plate at a temperature —T 0y on an arbitrary tem- 
perature scale. The direction of gravity is specified by 
the unit vector — k. The equations appropriate for our 
system are the Boussinesq approximations to the 
Navier-Stokes equations (Chandrasekhar, 1961, p. 16). 
We shall write these equations in a form in which the 
velocity and temperature (v and T) as well as the co- 
ordinate and time (r and /) are nondimensional. The only 
parameters entering the equations are then the Rayleigh 
number, R, and the Prandtl number, cr. The equations 
are 

Vv=0, (1) 

/I d \ 1 

V 2 ]V 2 v=-VX(VX(v-Vv))+^VX(VXkr), (2) 

\<r dt / <r 

d \ 

— v 2 )r=-v-(vr). (3) 

dt ) 

Equation (2) is actually the double curl of the mo- 
mentum equation, and hence the pressure variable is 
absent. The nondimensional variables are related to the 
dimensional ones (denoted by primes) in the following 
way: 


July 1963 


J. R. HERRING 


f 


327 



K 


T=T'/d , 
r=r'/d, 


K 

d 2 

Here k is the thermometric diffusivity of the fluid. 

The boundary conditions on the velocity field are 
derived from the requirement that the fluid exert no 
shear on the confining plates. This, together with the 
continuity equation, implies that all even derivatives of 
the vertical velocity, w, vanish on the boundary. In the 
nondimensional notation the boundary conditions are 

d m d m 

w(Q,t) = ?e>( 1,/) = 0, w = 0, 2, 4, * * • (4) 

dz m dz m 

and 

T(0, 0 = 0; 7X1,0= -1. (5) 


There are two more equations, for the x- and y-com- 
ponents of the velocity field, but these are not necessary 
for our problem. The last terms in the equations above 
for w and 8 have the form of a deviation of a bilinear 
fluctuating quantity from its horizontal mean (fluctuat- 
ing self-interaction). By discarding these terms we 
obtain the system to be investigated; 



where 


tf(s)=l 

dz 


d 

f. 

dz 


b) Discard of the fluctuating self -inter actions. It is con- 
venient to resolve the temperature field into a horizontal 
mean plus a fluctuating part; 

T— —z+\f(z } t)+6(Td). (6) 

Here, yp{zf) is a horizontally averaged distortion of the 
conduction state and 0(r 7 i) is the fluctuation of the 
temperature from its distorted value. In view of the 
boundary condition (5), and the interpretation of 8 as a 
fluctuation from the horizontal mean, we may write 

*(0, 0=*(1, /) = 0 (7) 

6(x,y,0,t) = 6(x,y, 1 ,/) = 0 = 0. (8) 

The bar on (8) indicates an average over the horizontal. 
We now introduce (6) into (1), (2), and (3) and subtract 
from each of the resulting equations their respective 
horizontal mean. We find 

n a \ 1 

( V* ) Vhv= RV 1 2 8+-{ vx vx (v • Vv) } „ 

\<t dt ) <r 

~~ V 2 V= ( (vd-kwd), 
dt ) \ dz) 

d d 2 \ d 

U— ^0, 

dt dz 2 ) dz 

where 

d 2 d 2 

Vi 2 = — + — . 

dx 2 dy 2 


The significance of omitting fluctuating self-interac- 
tion can be expressed formally by examining the 
hierarchy of moment equations obtained from (l)-(3). 
By multiplying equations (2) and (3) by u(0 T(t') 

and ensemble-averaging the appropriate sums of the 
resulting equations, we obtain the time evolution equa- 
tions for the correlation coefficients (i wf)> and (iuT'), 
and (XT'). These equations couple the above second 
order moments to the transfer terms, which are cubic 
in v and T. 

Since the system contains a non-vanishing first-order 
moment, the transfer terms contain both correlated 
third-order moments (cumulants) and products of first 
order moments with second-order moments. The dis- 
carding of the fluctuating self-interaction then corre- 
sponds to closing the system of moment equations by 
discarding the third order cumulants. 2 In the absence of 
mean fields this procedure would be empty. 

We must now verify that our procedure of deleting 
third-order cumulants does not lead to physically un- 
realistic results. For our procedure to be acceptable, the 
system of (9), (10) and (11) must obey the conservation 
laws associated with the complete set of convection 
equations, and they must be free from unphysical con- 
sequences of the sort recently discussed by Ogura (1962), 
for a similar problem in isotropic turbulence. Ogura has 
demonstrated that the assumption of zero fourth-order 
cumulants and nonzero third -order cumulants is in- 


2 Discarding third order cumulants is quite different from dis- 
carding third order moments. The latter procedure has as a conse- 
quence that no steady state nontrivial amplitudes exist. For an 
investigation of the dynamics of decay for zero third-order 
moments see Deissler (1962). 
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compatible with a positive energy spectrum for all 
wave numbers. 3 

With regard to the last point, it should be noted that 
the positive definite character of the kinetic energy wave 
number spectrum and the spectrum for the square of the 
temperature field follows directly from the fact that it is 
possible to write the equations which delete third order 
cumulants in terms of amplitudes rather than moments. 
We observe that the amplitude (9), (10) and (11) all 
have real coefficients; hence, the square of any ampli- 
tude will remain positive for all time if the amplitude is 
initially a real number. 

The conservation of entropy and kinetic energy are 
also preserved without the fluctuating self-interactions. 
By multiplying (10) by 0, (11) by \p, and adding, we 
obtain after integrating over the entire volume of the 
fluid, the equation of conservation of entropy, 

1 d 

{\ 0 \ + m»+{\W\ 2 +\W\ 2 }v={wd} v . ( 12 ) 

2 dt 

Here the v subscript indicates an integration over the 
entire volume of the system. 

We observe that equation (12), with a corresponding 
one for the conservation of the kinetic energy (Malkus 
and Veronis, 1958, p. 228) of the flow are exactly the 
same as those with the fluctuating self-interaction in- 
cluded. Contributions from the latter may be reduced 
to surface integrals which vanish. 

c) Fourier decomposition of the equations. It is con- 
venient to work with the Fourier components of (9)~-(ll) 
rather than their space-variable form. The free bound- 
ary conditions make the sine series appropriate. We 
therefore write: 

w(r,t)= Z f a (x,y)w n a sin«7rs, 


0(r,/) = Z M x ,y)0n a sin«7rz, 

n,m 

^(z,/) = Z sinmrs. 

n 

Here f Q (x,y) is an arbitrary set of orthonormal functions 
generated by the operator Vi 2 , and obeying appropriate 
periodic boundary conditions in the horizontal : 


/Id \ Aa 2 

+n 2 +a 2 )«,«=- OS, (13) 

\<r St / n 2 -\~a 2 

~+n 2 +a 2 je n a 

= ^n a Z #p(Wn+p a +<r(« — /O^^-pf), (14) 



7 rn 


— L E W p a {8 n+p*~\- v(p n)8 a \n-p\), (15) 

2 p=l a 


where 


and 


A=i0r 4 , 
r = 7 r 2 /, 

W n = WnA 2 , 

(t(oc) — 1 , #>0 

= 0, x~0 

= — 1, x<0. 


Manipulation of the convolution terms in (14) and (15) 
is aided by the following identities: 


Z A p(B„+ p -\-a(n p)B\ n _ p \) — Z^pC^N-pI ^4«4-p)? 

p p 


and 


Z A p (B n + p +<r(p n)B\ v „ n \) 

p 

= Z B p (An+p-\-<r(p— n)A jp-nj). 

p 

There are two conservation equations derivable from 

(14) and (15). The first is the Fourier representation of 
(12) for conservation of entropy. The other is the equa- 
tion that partitions the total heat flux between con- 
duction and convection; it is derived by multiplying 

(15) by \/n and summing over n. We find 

/Id \ 7T 2 

z +lk = -I«»V, (16) 

n \n 2 dt / 2 n.a 


Vi 2 /a(»,J')= -ir 2 a 2 f a (x,y) 

and 

Introducing the above representation into (9), (10) 
and (11) gives the following set of equations for the 
amplitudes w n , 8 n , and \p n : 


where 


Pn= — 7 T n ^n. 


Here the /Vs are the cosine transform of the mean tem- 
perature gradient. In the statistically steady state, (16) 
is the equation for the total heat flux, which is a con- 
stant of motion for the system. We now introduce a 
quantity N(t ), the total heat flux at the lower boundary: 


tf(O-i+Z0»(f). 

1 


3 For a complete discussion of the cumulant discard approxima- 

tions, see Kraichnan (1962a). 


( 17 ) 
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Fig. 1. Time development of the total heat flux, N(r) for 
R — 4X10 3 , 10 4 , 10 5 and a=1.5. The system is in the conduction 
state at t— 0, with all fluctuating amplitudes except W\ equal 
to zero. 


If the mean amplitudes are constant, 

7 T 2 

W=H £ Un a 8n a . 

2 n - a 

In our units, the conduction state transports unit heat 
flux and this equation is the nondimensional form for 
the familiar equation for the total heat flux. 

d. Structure of the equation . Before proceeding to the 
numerical results, we give a brief resume of the pertinent 
qualitative features of the system defined by (13), (14) 
and (15). First of all, we note that the horizontal wave 
numbers, a , are coupled only in their effect on the mean 
temperature field x//. This interaction occurs diagonally 
in the sense that each a interacts only with itself. As a 
consequence there is a degeneracy in the horizontal plan 
form of the motion; the system is insensitive to the 
particular cell shape. Moreover, the number of as is 
also indeterminate. The simplest situation is to have a 
single a support the motion and we investigate only 
this case here. 

A single a will give nontrivial answers for the ampli- 
tudes w and 8 only if it lies within a certain range. The 
range of a which will not support convection is obtained 
by assuming w and 8 to be small, and demanding that 
they subsequently decay. If w and 8 are small, \p will be 
small to second order and our question is equivalent to 
that of marginal stability (See Chandrasekhar, 1961, 
p. 35). The system then will not support convection if 

(1+a 2 ) 3 R 

>- (18) 

a 2 7 r 4 



Conversely, we assume that the steady state values of 
w and 8 will be nonzero if a lies in the range comple- 
mentary to (18). 

The time behavior of the system is complicated by 
nonlinear effects. In the approach to the steady state, 
our numerical results indicate that the system executes 
overdamped oscillations with an ever increasing period 
of oscillation. This last remark is understandable since 
w and 8 become marginally stable as t — > . 

If the mean field, \p, is statistically steady as t — » oc , 
we may use a theorem of Spiegel (1962, p. 196) to show 
that w and 8 are independent of time. Spiegel has shown 
that the principle of exchange of stability is valid (for 
free boundaries) in the presence of the mean gradient 
corresponding to the steady state solution to the mean 
temperature field given by (15). This implies that the 
growth rates for the appropriate eigen-function expan- 
sion for w and 8 must all be zero in order for there to be 
a statistically steady state. 

3. Numerical procedure 

In performing the numerical integration, we discard 
from the onset those Fourier amplitudes which will be 
zero in the steady state. We assume that the steady 
state amplitudes to, 9 and 0 have even parity about the 
mid-point z—\. This means that the even sine modes of 
co and 8> and the odd cosine modes of /3 will have zero 
amplitude in the steady state. We therefore put their 
initial values equal to zero. The equations of motion 
(13), (14) and (15) then imply that the odd parity modes 
will remain zero for all subsequent time. Defining 
j3n=ftn, we may rewrite (14) and (15) in a more con- 
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Fig. 3(a). Mean temperature, T(z ), for R = 4X10 3 and <* = 0.8. 

Fig. 3(b). Mean gradient, / 9(z), for R = 4X10 3 and c* = 0.8. 0(z) is normalized by the total heat transport, A = 3.92. 



venient form: 



— 0?«+^ 2/>)cO]„_2pj), (14') 

p-i 

d \ _ 

— +4w 2 )/5 n 
dr / 

= — 2tT 2 W 2 £ C0 J ,(^2/ l -p+O'(/)~2w)^|2n-p!). (15') 

P“1 


Here, we have dropped the a superscript since we are 
interested in the system containing only one a. Equa- 
tion (13) remains unchanged and the total heat flux is 
computed from (17). 

Our procedure for integrating these equations is to 
assign an initial set of amplitudes to w n , 6 n , /3„, and . 
allow the system to evolve to the steady state. In doing 
so, we must truncate the infinite set of equations. Our . 
procedure in this matter is to set all amplitudes co n , 0 ny 
fin of index greater than a certain integer, n Q , equal to 
zero. This method of truncation guarantees exact con- 
servation of heat flux and entropy for the abbreviated 
system. Since an is generally large and the 0„’s decrease 
rather slowly, we see from (15') that amplitudes for fi n 
above no/ 2 will have significant truncation errors. 
Truncation errors are assumed to be negligible if in- 
creasing no does not appreciably alter the value of the 
total heat transport. The total number of p N modes 
included in these calculations ranged from 20 modes at 
R = 4000 to 80 modes at R= 10 6 . The errors in the total 
heat transport due to the above are estimated to be less 
than one part in 10 3 . 

The integration forward in time was continued until 
constancy of heat flux (16) and entropy (12) was 
achieved to one part in 10 4 . The time, in r units, neces- 
sary to achieve this ran from ~ 1.4 at R = 4X10 3 to 0.3 
at R= 10 6 . At high Rayleigh numbers, this criterion was 
not too satisfactory, since constancy of heat flux and 
entropy were achieved long before the amplitudes w 
and 6 became steady. For these cases, it was necessary 
to check the time derivatives of the slowest evolving 
amplitudes, w\ and The system was observed to be 
steady if the derivative of w\ was less than 2 per cent 
of W\. 
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Fig. 5(a). Mean temperature, T(z), for R= 10 4 and a= 1.0. 

Fig. 5(b). Mean gradient, jS(z), for R=10 4 and a=1.0. (3(z) is normalized by the total heat transport, #=5.82. 


Examples of the time evolution for the total heat 
transport N(t) are given in Fig. 1, for R=4X10 3 , 10 4 , 
10 5 and a= 1.5. The system was started in the conduc- 
tion state at r = 0, with all fluctuating amplitudes co n 
and On equal to zero, except c*u, which has an initial 
value of unity. The convection is seen to develop initially 
by way of large oscillations, and to decay to the steady 
state with overdamped oscillations, whose period be- 
comes increasingly larger. The time scale of the initial 
oscillations in these curves is of the order of the growth 
rate time in the conduction state. 

4. Discussion of results 

The computed steady state amplitudes are shown in 
Figs. 2—13. The normalization for 0, w and 6 is given in 
the captions, while T{z) requires no normalization. The 
graphs of T{z) are in a reflected coordinate system to 
conform with an accepted procedure. The values of a in 
Figs. 2-9 were chosen so that the heat transport is very 
near its maximum. We now discuss in some detail the 
physical features of the steady state amplitudes co, 0, 
0 and T. 

a) The mean temperature T , and mean gradient 0(z). 
The mean fields, T and 0 in Figs. 3, 5, 7, 9, 11 and 13 
have an interesting behavior near the boundaries. At 
low Rayleigh number, these fields closely resemble the 
perturbation results of Markus and Veronis (1958), but 
the temperature gradient is slightly negative in the 
central region. As the Rayleigh number is increased, the 
negative temperature region collects closer to the 
boundary while in the central region, the temperature 
gradient becomes extremely small but positive. 



The negative temperature gradient boundary region 
is apparently produced by an overshoot phenomenon. 
These occur typically for motions of large horizontal 
scale (small a) and disappear for motion of small hori- 
zontal scale. (See Figs. 9, 11 and 13.) If the motion has 
a large horizontal scale, an element of fluid close to the 
lower plate moves in a region of high temperature for a 
relatively long time. When it eventually turns upward, it 
moves unchecked by eddy processes and penetrates the 
body of the fluid with an excessive heat flux. The mean 
gradient accommodates this motion by turning negative. 
The negative 0 region then checks the velocity field, so 
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Fig. 7(a). Mean temperature, f{z ), for R=10 6 and a =1.5. 


Fig. 7(b). Mean gradient, /3(z), for R= 10 6 and a = 1.5. /S (z) is normalized by the total heat transport, N= 13.82. 



that the advective heat transport decreases toward the 
middle of the fluid. We see here evidence for the non- 
local property of the flow; if the Rayleigh criteria for 
convection were applicable locally, a negative /3 region 
would not persist in the steady state. 

For motions of small horizontal scale (Figs. 11 and 13) 
the situation is somewhat different. In this case, an 
element absorbs little heat from the lower boundary 
region and loses it quickly by conduction because it 
belongs to a vertically elongated cell pattern. It also 
loses momentum by viscous drag, and attains its 


terminal velocity before reaching the central region of 
the fluid (see Figs. 10 and 12). To maintain constancy 
of heat flux the central region must conduct rather 
strongly, so that the mean gradient becomes large there. 

b) Velocity and temperature fluctuations. The velocity 
and temperature fluctuation fields are shown in Figs. 2, 
4, 6, 8, 10 and 12. We observe that the velocity fields, 
for all Rayleigh numbers, have an extremely large first 
mode. For example, at R = 4000 (Fig. 2) w x represents 
99 per cent of the total velocity amplitude, while at 
R= 10 6 (Fig. 8) Wi is 95 per cent of the total. On the 
other hand, the B n modes decrease rather slowly as 
n increases. 

The above behavior of the co n and 6 n spectra displays 
the character of the nonlinear coupling in our system. 
Thus, the term tends to be the dominant con- 

tributor to p n [see equation (15')] for reasonably small 
n. Conversely, terms of the form jS n coi and tend 

to be the dominant contributors to 6 n [equation (14')]. 
The nonlinear coupling scheme in the equations of 
motion is therefore highly nondiagonal, as opposed to 
the case of isotropic turbulence. 

The strong nondiagonal coupling in the system of 
Fourier modes is a result of the distortion of the mean 
temperature profile combined with the pressure and 
dissipative forces for incompressible flow. The above 
forces are directly responsible for the occurrence of 
sixth -order derivatives in the marginal stability prob- 
lem, of which the steady state amplitudes co and 0 
are solutions in the presence of the mean field (3. If we 
solve for the velocity modes c in the presence of the 
mean gradient /3, by using the iteration technique of 
Section 5, we see that the higher modes of u n are sup- 
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Fig. 9(a). Mean temperature, ?(z), for R— 10 6 and <*=1.5. 


Fig. 9(b). Mean gradient, /3(z), for R= 10® and <*=1.5. /3(z) is normalized by the total heat transport, iV=31.48. 


pressed by a factor ^n~ B . For a reasonable 0, this factor 
results in the higher co„ modes making only a small 
contribution to a>. 

c) Temperature gradient spectrum . The cosine spec- 
trum of the mean temperature gradient, 0 ,is given in 
Fig. 14 for R= 10 4 , 10 5 , 10 6 and a=1.5. We have con- 
. nected the points with a smooth curve for the sake of 
clarity. We notice a tendency for the lower modes to 
saturate at 0»=2, which corresponds to the small 
gradient outside the boundary layer. In fact, if fi n —2 
for all n , fi(z) is a 8 function, and the gradient vanishes 
everywhere except at the boundary, where it becomes 
singular. At large Rayleigh numbers, the fi n spectrum 
is nearly Gaussian for small n , but decreases more 
rapidly at large n. 

The tendency for the ft’s (for small n) to approach 2 
as an upper bound is closely connected with the fact that 
the velocity field is marginally stable on the mean 
temperature gradient, 13. This feature is brought out 
more clearly by examining the relation connecting the 
mean gradient field, ft and the Rayleigh number R. 
Using the iteration method of Section V, we find 


(1 + 2 ) 3 7 T 4 

a 2 R 




1/4 


(1+a 2 ) 3 


1-^ft i [(2*+l) 2 +a 2 ] 3 


i ' 1 r 



I i i i i 

0.2 0.4 0.6 0.8 Z 

Fig. 10. 3.22X10“ 3 w> and 8.570 for R = 10 6 and <* = 6.0. 


The computed spectra (Fig. 14) are qualitatively 
quite different from the one derived by Malkus (1954, 
p. 200). His spectrum is given by 


X(0n-0n +l ) 2 +**‘. (19) 

This series for R~ ! converges rather rapidly for all the 
ft which have been computed, and the terms explicitly 
written in (19) give R to an accuracy of < 20 per cent at 
R=10 6 . We note that for this equation to balance at 
large R, ft must approach 2, and the remaining lower 
modes must decrease rather slowly as n increases. 


ft^fl — Y 

\ 2n 0 +lJ 

Here 2w 0 +l is a cut-off in the fi n spectrum, and it is the 
total heat flux in our units. Marginal stability is 
achieved at a much lower Rayleigh number for this 
spectrum than for the ones computed here. 
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Fig. 11(a). Mean temperature, f(z ), for R=10 6 and oc-6.0. 

Fig. 11(b). Mean gradient, /S(; z), for R = 10 6 and a = 6.0. /S(z) is normalized by the total heat transport, iV = 22.3. 



With regard to the Malkus theory, Figs. 10 and 11 are 
relevant. For this case (R= 10 6 , a = 6.0) the temperature 
gradient is everywhere positive except near the bound- 
aries where it approaches zero. The fields in Figs. 10 and 
11 therefore fulfill all the requirements of the Malkus 
theory as formulated by Spiegel (1962). We note for 
this case that the total heat flux is ^22, whereas Malkus 
obtains a heat transport of ~11 for free boundaries. In 
making this comparison, one should remember that 
these computations were made for a single horizontal 
wave number, whereas Malkus presumably allowed 
for a full spectrum of as. However, if we interpret the 
computed heat transport as an upper bound to the 
heat transport as the Malkus theory prescribes, we con- 


clude that for free boundaries the actual upper bound 
is at least a factor of two larger than that obtained by 
Malkus. 

d) The total heat transport as a function of a and R. 
The total heat transport, as a function of R and a is 
given in Fig. 15. The Rayleigh numbers are indicated 
in the figure. These curves closely resemble the pertur- 
bation calculations at small R, but become increasingly 
broadened as the Rayleigh number increases. For a 
given R, the heat transport is entirely conductive 
(N= 1) if a lies outside the bounds prescribed by (18). 
The value of a which maximizes the heat transport is 
1/V2 at the critical Rayleigh number (R=657), and 
apparently increases linearly in R*, at large R. It is well 
represented at high Rayleigh numbers by the formula: 

0.7+0.01R*. (20) 

The data on this point is not entirely conclusive because 
of the large breadth of the curves. It should be pointed 
out that (20) cannot be a correct asymptotic formula 
since a m ax is proportional to R* and the value of a 
beyond which a single a cannot support convection is 
proportional to R* (equation 18) . 4 An estimate of the 
Rayleigh number beyond which (20) is incorrect is not 
warranted by the accuracy of the curves, but according 
to (18) and (20), it isR>10 18 . 

The maximum heat transport as a function of R is 
given in Fig. 16. For R>3000 the data are accurately 
represented by the following R* law: 

iV=0.31Rk (21) 


4 R. H. Kraichnan has pointed out to the writer that the heat 
transport, A max , is then asymptotically proportional to R 3/10 . 
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Fig. 13(a). Mean temperature, T (z), for R=10 5 6 and a = 9.0. 

Fig. 13(b). Mean gradient, 0 ( 3 ), for R= 10 6 and a = 9.0. £( 2 ) is normalized by the total heat transport, N= 5.40. 


Experimentally, the Nusselt number N is ^ 0.085 R*, 
for large R, and rigid boundaries (Jakob, 1959). We see 
no evidence for an intermediate R* law, but such a law 
may only be obtained in the rigid boundary problem. 
Below R^IO 3 , the data fits smoothly to the perturba- 
tion calculation of Malkus and Veronis (1958). 

The discrepancy between (21) and the experiment is 
partly a result of eddy processes which our procedure 
omits. The use of rigid boundary conditions may im- 
prove the agreement, but if (21) is corrected for 
boundary effects as done by Malkus (1960) by de- 
creasing N by (657/1 107)*, there remains a discrepancy 
of a factor of 2. If we choose horizontal wave numbers 
such that the mean gradient is everywhere positive 
(Figs. 10 and 11) the discrepancy is reduced to 1.8. The 
latter fields, however, have the unattractive feature of 
having a large temperature gradient in the central 
region of the fluid. 

5. Linear stability of the fields 

The velocity and temperature fields we have so far 
discussed are stable against the introduction of a dis- 
turbance of the same horizontal wave number for which 
the fields were computed. This stability is inherent in 
the method of integrating the equations. The stability 
of the steady state amplitudes against disturbances at 
wave numbers a other than that a which supports the 
convection process has not yet been assured in our 
calculations. The question of stability of the solutions 
against disturbances of finite amplitudes leads directly 
back to the multi-a system of equations (13), (14) and 
(15). We should assume a whole spectrum of a’s are 
initially excited, let them evolve to steady state, and 



Fig. 14. Cosine spectrum of the mean temperature gradient 
for R — 10 4 , 10 6 , 10 6 and a= 1.5. 


repeat the calculation for an ensemble of initial condi- 
tions. We shall be content here with an investigation of 
linear stability of the system. This problem has some 
intrinsic interest, but our main purpose is to lay the 
framework for an investigation of the multi-a system. 

It is convenient to pose the linear stability problem 
in terms of the Fourier amplitudes [equations (13), (14) 
and (15) J. We suppose that the system w n a \ 6 n ai and 
p n have their steady state values, introduce disturbances 
8co n a , 8d n a and 8/3 n and ask whether the latter grow or 
decay. Since a 1 and a are not coupled, must decay 
initially. The problem then is reduced to determining 
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Fig. 15. The total heat transport N as a function of a for 
R = 4X 10 3 , 10 4 , IQS 5x105 and 10 6 . 


the growth rates for 8w n a and 86 n a in the presence of the 
mean gradient ft n . Since exchange of stabilities has been 
proved for this system (Spiegel, 1962) we know that the 
system w n ai , 6 n ai , and ft n will be stable if the smallest 
critical Rayleigh number R c for the perturbation sys- 
tem, 8oo n a , 86 n a is larger than the Rayleigh number for 
which w al and 6 al were computed. 

The marginally stable amplitude 8w n a , and 80 n a 
satisfy (13) and (15) at a Rayleigh number R c , with the 
time derivatives put equal to zero. Since the smallest R c 
takes an eigen function even about z—\, we may ab- 
breviate the perturbation system by eliminating the 
even sine modes from the velocity and temperature 
fluctuations. Defining 

<p n ^ 8W2n-h 
< Pn ~ 8W 




Fig. 17. Critical Rayleigh number R c for R = 5X10 6 as a func- 
tion a. The value of a which supports the mean temperature fields 
labels the various curves. 


we may eliminate 86 n by using the steady state form of 
(13) and write the marginally stability problem in the 
following matrix form : 

A(P. )<p=h<p, (22) 

where 

a 2 

A (ft) = ~~ — -{ 8nm+i(ft\ n—m\ ftn+m- 1 )}, 

[_(2n-iy+oi 2 J 

7T 4 

jU= — , /3o=0. 


In writing the matrix A , we have used the alternative 
form for the convolution term in (14). 

The largest eigen value, m (smallest R c ), of (22) may 
be obtained by the matrix iteration technique (Hilde- 
brand, 1952). Since the first sine mode of the velocity 
will be largest, we may conveniently begin the iteration 
on a vector containing only this mode. Defining 


1 1)= (1.0,0,- • -,o,- • •) 

we may write 

7T 4 _ (lM n, ,l) 

~ jUrnux ~ lim . 

R c ( a) <1 M”- 1 ! 1> 


(23) 


The convergence of the iteration scheme is quite rapid 
because of the structure of the A matrix. At the highest 
Rayleigh number considered, R= 10 6 , the 11th iteration 
gives R c to one part in 10 6 . 

The calculated R c (a) , s are shown in Fig. 17 for 
R=5X10 5 . The particular a which supports the mean 
field, ft, labels the various curves. At low Rayleigh 
numbers (R<5.10 3 ) these curves closely resemble those 
produced by perturbation calculations. Above R^IO 5 
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they become increasingly distorted; steady state ampli- 
tudes of small horizontal wave numbers are enormously 
unstable with respect to an introduction of a disturbance 
at large a. 

In Figs. 18 and 19 we give the zones of instability for 
the computed amplitudes for R=10 4 and R=5X10 5 . 
In these graphs, a\ is the wave number that supports 
the convective process, and a 2 is the wave number of 
the perturbation amplitudes. The regions of instability 
are indicated by the shaded areas, whose outer bound- 
aries are lines of marginal stability. The line ofi = 0:2 is a 
trivial case of marginal stability. The value of a at which 
the two curves cross represents a solution which is 
infinitesimally stable against all other as. This value 
of a begins at 1 /V2 at the critical Rayleigh number and 
increases slowly with increasing R. The rate of increase 
is seen to be slower than that a which maximizes the 
total heat transport. Referring to Fig. 15 we see that 
the use of the most stable a instead of a ma x will not 
appreciably change the total heat transport. 

The zones which linear stability theory predicts must 
have two or more a } s supporting convection are indi- 
cated by the cross-hatched regions shown in Figs. 18 and 
19. These regions are obtained by perturbing the a 1 
fields at 0 : 2 , assuming that the <*2 field subsequently 
dominates the convection, and then demanding that the 
a* fields be unstable with respect to a perturbation at ol\. 
The cross-hatched region is then bounded by the de- 
' scending marginally stable curve and its reflection about 
the 45 deg line. At small Rayleigh number, R<4X10 3 , 
• this area vanishes but it gradually increases with 
Rayleigh number. 

6. Concluding remarks 

The temperature and velocity fields computed here 
with the fluctuating self-interactions absent show quali- 
tatively a reasonable behavior. The boundary-layering 
of the temperature field, which is found experimentally, 
is faithfully reproduced by the system, and the heat 
transport has the experimentally determined dependence 
on Rayleigh number. In this respect, our results for the 
velocity and temperature amplitudes, as well as the 
stability analysis of the fields, confirm the original ideas 
of Malkus. However, our result for the heat transport 
for free boundary conditions does not agree quantita- 
tively with Malkus. 

The only disquieting features of the results are the 
negative temperature gradients which can occur near 
the boundary for small a, and the rather large amount 
of heat transported by the system. Aside from eddy 
processes, there are two other modifications in the 
system which must be explored before its quantitative 
accuracy can be properly accessed. 

First, the use of the more realistic rigid boundary 
conditions will enable one to examine quantitatively the 



Fig. 18. Stability diagram for R= 10 4 . ai is the wave number 
that supports convection, and <22 is the wave number at which a 
small perturbation is introduced. The shaded region indicates 
instability. 



Fig. 19. Stability diagram for R = 5X 10 6 . ai is the wave number 
that supports convection, and a 2 is the wave number at which a 
small perturbation is introduced. The shaded regions indicate 
instability. 

role of the eddy processes in producing the experimental 
temperature profile and the total heat flux. The presence 
of shear forces at the boundary will decrease the com- 
puted heat flux, and in checking the development of 
large scale horizontal motions there, it will reduce the 
negative temperature gradient. Preliminary indications 
are that the use of rigid -boundary conditions decrease 
the total heat transport by a factor of 2.3. Secondly, the 
introduction of several horizontal wave numbers will 
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make the system more realistic, particularly at large 
Rayleigh numbers. It will also permit a study of finite 
amplitude stability of the system. The above modifica- 
tions are currently under investigation and will be 
reported on in the near future. 
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